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I. INTRODUCTION 



Over the last two decades, discrete Euclidean path integrals constructed using the so called Trotter decomposition!]]'! 
have been widely used as the starting point foe quantum Monte Carlo (QMC) simulations of lattice models at finite 
temperatureJ3i9 In these "worldline" methodsjj the discretization At in imaginary time introduces a systematic error 
in computed quantities, which in principle can be eliminated by carrying out simulations for different discretizations 
and extrapolating to At = 0. The error typicallyem scales as (At) 2 . However, algorithms can also be constructed 
which are associated with no inherent systematic errors, thus eliminating the need for multiple simulations and 
extrapolations. Such a QMC scheme, applicable to the ferromagnetic spin-1/2 Heisenberg model, was first devised 
Ch ' by Handscomb as early as in 1961E This method is based on a series expansion of the density matrix operator 
exp(—(3H), which in the case of the Heisenberg model can be written in terms of products of permutation operators. 
i Their traces can be evaluated exactly and are positive definite. One can then carry out importance sampling in the 

space of operator sequences, and obtain results exact to within statistical errors. Handscomb's method is not directly 
t-H , applicable to other models Jj not even for the antiferromagnetic Heisenberg model for which the scheme breaks down 
due to the non-positive-definiteness of the traces.Ej There was therefore not much follow-up on Handscomb's pioneering 
efforts, and there was li±±le progress towards practical QMC algorithms until Suzuki proposed the use of the Trotip: 
formula in this context.H Several methods based on this controlled approximation were subsequently developedfoEJ 
Variants of Handscomb's method were also later developed, for the antiferromagnetic Heisenberg model by Lee et 
aZ.Jl3 and for the XY-model by Chakravarty and Stein,ll3 but this type of simulation scheme was fac,.a_.long time 
£ — ' still perceived as fundamentally limited by its reliance on the special properties of spin-1/2 operators .EjEj However, 
Q\ . the generalizations of Handscomb's method developed for S > 1/2 spin models by Sandvik and Kurkijarvi,E3 and 
for ID Hubbard-type models by Sandvik,U have now clearly demonstrated that non-approximate algorithms based 
on "stochastic series expansion" (SSE) can in principle be constructed for any lattice model. In this scheme, the 
basic limitation of the earlier formulations of Handscomb's method is overcome by expanding the traces as diagonal 
L | matrix elements in a suitable chosen basis. The importance sampling is then carried out in a space of basis states 
and operator sequences. E£(liH In practice, this type of method is of course still limited to models for which a positive 
definite weight function can be achieved. The cases for which this is possible coincide with those for which the weight 
is positive definite also in the standard Trotter-based path integral formulations. In fact, despite the different starting 
points of the two approaches, the SSE configuration space is strongly related to an Euclidean path integral. Many 
. , characteristics are therefore shared, including the range of practical applicability, the types of observables accessible 
for evaluation, and the scaling of the computation time with the system size and the inverse temperature (3. The 
5-H \ main advantage of SSE is of course the absence of systematic errors. It should be noted that in order to eliminate the 
Trotter error in worldline calculations, simulations iave to be carried out for several sufficiently small discretizations 



At. Since the computation time scales as 1/(At) 3 ,L3 the SSE approach can in practice be considerably more efficient 
than the worldline method in cases where completely unbiased results are needed. 

Recently, other approaches to exact QMC algorithms have been proposed. Beard and Wiese. succeeded in formu- 
lating a worldline algorithm for the spin-1/2 Heisenberg model directly in the At — > limit £3 Constructed within 
the framework of a non-Metropolis sampling scheme with global "loop-cluster" updates previously developed by Ev- 
ertz et aZ.o (a generalization of the classical cluster spin algorithmEl), this method also has the added advantage of 
significantly shorter autocorrelation times. Concurrently, Prokof'ev, Spristunov and Tupitsyn suggested the use of the 
standard perturbation expansion as a basis for a QMC path integral.tSI For a finite system at finite temperature the 
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series converges to an exact result for a finite number of terms. A scheme involving a novel class of local updates was 
suggested for efficiently sampling the continuous-time paths. ■— ■ 

It is clear that the methods by Beard and Wiesalia and Prokof'ev et al&3 are strongly related to each other, involving 
the same configuration space of world lines in continuous imaginary time, but differing in the , sampling procedures. 
Continuous-time path integrals also have many properties in common with the SSE sum.l 16 H 17 l Notably, a transition 
event in imaginary time corresponds directly to the presence of an off-diagonal operator in the SSE operator string. 
Here we discuss this connection in detail, and introduce a simple modification of the SSE algorithm for simulations 
in the interaction representation. This formulation can be expected to be more efficient than standard SSE in cases 
where the diagonal part of the Hamiltonian dominates. In order to explore the properties of the new method, we study 
the spin-1/2 Heisenberg chain, as well as a spin chain including couplings to dynamic (fully quantum mechanical) 
phonons. We consider a coupling via a linear modulation of the spin exchange by a local dispersionless oscillator 
(Einstein phonon). The system undergoes a spin-Peierls (dimerization) transition at zero temperature. A detailed 
study of the model, and its relevance for understanding the magnetic properties of the recently discovered spin-Peierls 
compounds GeCuC>3 (Ref. |||) and a'— NaX^Os (Ref. ^3|), will be presented elsewhere.cil Here we only consider a 
single set of model parameters, in the general regime expected to be of physical relevance, and illustrate the use of the 
new method to calculate a variety of physical observables, both at finite temperature and in the limit T — > 0. Based 
on the results, we conclude that the effects of dynamic phonons cannot be neglected in quantitative descriptions of 
materials such as those mentioned above. 

The outline of the rest of the paper is the following: In Sec. II we review the formalism of the SSE method. The 
perturbation expansion in the interaction representation and its relation to the SSE scries are discussed in Sec. III. In 
Sec. IV we implement an interaction representation algorithm for the Heisenberg chain, and discuss the performance 
of the method. In Sec. V we consider the spin-phonon model. Readers interested mainly in the new results for this 
model are advised to skip directly to the introductory part of Sec. V, and then go directly to the results section V-B. 
The discussion there does not rely heavily on the previous, more technical parts of the paper. Sec. VI concludes with 
a summary and outlook for further developments and applications of the new QMC algorithm, in particular to various 
models including phonons. 



II. STOCHASTIC SERIES EXPANSION 



Here we review the general formalism of the SSE method, needed as-,a basis for the discussion in the following 
sections. More details of the algorithm have been described elsewhere.E3 Some recent applications to spin systems 
and ID fermion systems are listed in Refs. |2^-|30|. 

The starting point for evaluating an operator expectation value at inverse temperature /3, 

(A) = |Tr{ie-^}, Z = Tr{ e -^}, (I) 

is to Taylor expand exp(— (3H) and write the traces as sums over diagonal matrix elements in a basis {|o?)}. The 
partition function is then 



("/?) 

a n=0 

The Hamiltonian is next written as 

M 



z = EE^^"i«)- ( 2 ) 



H = ^H b , (3) 



b=l 

where the operators H b have the "non-branching" property, 

H b \a) =h b {a,P)\P), (4) 

where |a) and are both basis states in the chosen representation. Each power H n is now expanded as a sum over 
all possible products of n of the operators H b . With S n denoting an index sequence referring to the operators in the 
product (the operator string), 

Sn = (&i,..., b n ), he {!,..., M}, (5) 
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the partition function becomes 



a n=0 S„ i=l 

Eqs. (||) and ([!]) of course represent a completely general formal device. In practice, one typically chooses the basis 
so as to make the sum (|J) as simple as possible. For example, for the Heisenberg model, 

ff = J^S i -S i , (7) 

the eigenstates of Sf can be chosen; \a) = \Sf, . . . , S^). Writing the Hamiltonian as 

H = J + ±(S+S- + Sr S +)] , (8) 

all the two-spin operators satisfy the requirement ([I]). For S — 1/2, S^SJ~ + satisfy (j|) and can be considered 

as a single operator, whereas for S > 1/2 the two terms have to be treated as separate operators. For tight-binding 
fermion or boson models, the real-space occupation number basis is typically chosen. For fermions and hard-core 
bosons the hopping operator cfcj + cttk satisfies (^), whereas for unconstrained bosons the terms again qualify only 
individually. 

For a finite system at finite fj, the lengths n of the operator strings contributing significantly to the partition function 
are restricted to a finite range. In a Monte Carlo simulation of the series expansion, terms (a, S n ) arc sampled with 
a probability proportional to the weight function corresponding to Eq. ([|): 

W(a,S n ) = ^^(a\f[H hi \a). (9) 
t=i 

Here it will be assumed that W(a,S n ) is positive definite, which of course is not.-always the case. With a non- 
positive-definite weight, simulations can in principle still be carried out using |VF|,B but in practice the statistical 
fluctuations of calculated expectation values diverge if the positive and negative contributions almost cancel each 
other, which-they do exponentially both with increasing system size and decreasing temperature (the infamous sign 
problem) £M3 As already discussed in the Introduction, this is the most severe limitation of the method— — shared 
also by standard techniques such as the worldline method (in the case of fermions, "determinant methods"liil typically 
are more effective in dealing with the sign problerrH). Still, the class of models for which a positive definite W can 
be achieved is significant enough to motivate the continuing development of more efficient QMC methods for their 
study. 

Proceeding as in the derivation of (||) , the numerator in ([!]) corresponding to a given operator A of interest is also 
expanded. If the expectation value can then be cast into the form 

{) E a E n E Sn W(a,S n ) ' (10) 

the simulation estimate of (A) is given by the average of the estimator A(a, S n ) over the sampled configurations: 

(A) = (A(a,S n )). (11) 

The formally simplest observable within the framework of the SSE method-is. the internal energy, E = (H). As in 
Handscomb's original formulation, the estimator involves only the power n;cHl3 

E = -(n)/0. (12) 

This, in combination with the expression for the heat capacity! 

C= (n 2 ) - (n) 2 - (n), (13) 

shows that the terms contributing significantly are of length ~ f3N (at low temperatures), where N is the system size. 
A derivation of (FL2|) will be discussed in the next section. 
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Already at the level of Eq. (||), the close relationship between SSE and a standard Euclidean path integral is evident. 
The operator string defines a set of propagated states \a(p)), p = 0, . . . , n: 

v 

\a(p)) ~l[H bi \a), \a(0)) = \a). (14) 
1=1 

A nonzero weight @ implies the periodicity condition |a(0)) = \a(n)). The propagation index p plays a role analogous 
to imaginary time in a standard path integral. The exact relation to imaginary time can be obtained by deriving 
an expression for a time-dependent correlation functionJ13 Consider two diagonal operators Ai and Aj. In a given 
configuration (a,S n ), their eigenvalues in the states are denoted a.i[p] and One can show that the 

time-dependent correlation function d 3 (t) = (Aj(T)Ai(O)) , where e rH Aje~ rH , is given bytH 



r - ir!= ^(l) T '" ( V r ^ w) - 



E 

*m=0 

Here (m) is a correlator between states separated by m propagations: 



1 " 

Cij(m) = — -^a,j\p + m]ai\p]. (16) 
n + 1 p=0 

The periodicity of the propagated states of course implies that ay [p + n] = cij [p] . In the equal-time case, only m = 
contributes to (|l5|), and CV, (0) is simply given by Oj[p]oj[p] averaged over p: 

C ij {Q) = (— 7l Y. a M^m)- (17) 

Eq. ( |l5| ) shows that an imaginary time separation r corresponds to a distribution of separations m between propagated 
states, peaked around m = tit/ (3. Hence, the propagation index p in the SSE method indeed is closely related to the 
time in an Euclidean path integral. 

As already discussed above, the range of contributing powers n is limited in practice. One can therefore explicitly 
truncate the series expansion at some maximum power n — L, large enough to introduce only an exponentially small, 
completely negligible error. By inserting L — n unit operators in the operator strings, a configuration space is obtained 
for which the sequence length formally is fixed at L. Defining the unit operator Hq = I, the summation over n in (|^) 
is implicitly included in the summation over all sequences Sx, if the range of allowed indices is extended to include 
also bi = 0. The weight (|^) has to be divided by (^), in order to compensate for the number of different ways of 
inserting the unit operators, resulting inllj 

W(a, S L ), { ~ m ^~ n)l (a\f[H bi \a), (18) 

i=i 

where n now denotes the number of non-0 indices in Sl- This fixed-length formulation is useful for the construction 
of an efficient sampling scheme for the sequences. For purposes of measuring operator expectation values, one can 
still use the expressions discussed above, with the sequences S n obtained by omitting all the zeros in the generated 
Sl- 

In order to ensure a sufficiently high truncation L, the power n is monitored during the equilibration part of the 
simulation. If n exceeds some threshold value L — A^, the sequence is augmented with, e.g., 2A^ randomly positioned 
unit operators, corresponding to L — > L + 2A^. With A^ i=s L/10, this procedure typically converges rapidly to a 
proper L. During a subsequent simulation (of practical duration), n never reaches L. The truncation is therefore no 
approximation in practice. 

The details of the Monte Carlo sampling procedures of course depend on the model under consideration. Here 
only some general principles will be discussed. The operators H a can be divided into two classes; diagonal and 
off-diagonal. There are no a priori constraints on the number of diagonal operators that can appear in Sl- The 
probability of a diagonal operator -f/dia at a position p is only determined by the state \a(p — 1)) on which it operates. 
The general strategy for inserting and removing diagonal operators is to attempt substitutions with the unit operator 
Ho introduced in the fixed-length scheme (note again that Hq is not part of the Hamiltonian) : 

H «-> H dia . (19) 
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This can be attempted consecutively at all positions in Sl- The weight change needed for calculating the Metropolis 
or heat-bath acceptance probability involves only the matrix element (a(p — l)|i?dia|a(p — 1)) and the prefactor 
(—(3) n (L — n)\, with n changing by ±1. With |a(0)) stored initially, the subsequent states can be generated one-by- 
one as needed during the updating process. 

Suitable constants have to be added to the diagonal operators in order to make all the eigenvalues of —f3Hdi& 
positive. According to Eq. (|l8|), the presence or absence of a sign problem then depends only on the off-diagonal 
operators H Q g . They are associated with various constraints, and cannot be inserted or removed at a single position 
only. They can always be inserted and removed pairwise. One way to do this is in substitutions with diagonal 
operators, according to 

Hdin, -ffdia «-> H g,H oS . (20) 

In some one-dimensional models, the above types of updates are sufficient for achieving ergodicity. In other cases, 
more complicated updates are also required (e.g., involving off-diagonal operators forming loops around plaquettes in 
2D). The constraints and weight changes associated with local updates involve only operators present in Sl which 
act on a small number of lattice sites surrounding those directly affected by the update. Typically, this allows for a 
sampling scheme for which the computation time scales as N[3x-H 



III. RELATION TO THE PERTURBATION EXPANSION 



In this Section we discuss the general principles of carrying out importance sampling of the standard perturbation 
expanaipn in the interaction representation. This starting point for a QMC scheme was recently suggested by Prokof'ev 
et al.x3 We here show that the configuration space of this method is closely related to that of the SSE method. We 
also derive expressions for several types of observables. 

The partition function for a Hamiltonian 

H = D + V, (21) 

with a diagonal (unperturbed) part D and an off diagonal (perturbing) part V is given by the standard time-ordered 
perturbation expansion in V , 

Z = J2(-l) n dn dr 2 --- dT n Tr{V(n)V{T 2 )---V(r n )}, (22) 

„=o J ° J ° Jo 

where the time dependence in the interaction representation is V(t) = e TD Ve~ rD . In the same way as was done in 
the SSE scheme, V can be decomposed into operators that satisfy the requirement (Q), now in the basis {\a}} where 
D is diagonal: 

M V 

V = J2H b . (23) 

6=1 

For a given model, the operators in the above sum are of course a subset of those in the SSE Hamiltonian ([?]), where 
we now define the indexing such that all with b > My are diagonal. An index sequence defining a product of n 
of the operators Hb is defined as before. In order to distinguish the SSE sequence S n , which contains off diagonal 
as well as diagonal operators, from the perturbation expansion sequence containing only off-diagonal operators, we 
denote the latter by T n : 

T n = (&i,..., b n ), b p e{l,...,M v }. (24) 
Expanding the trace over diagonal matrix elements gives 

Z = E £ E f dTl / 1 dT2 ■ ' ■ / " 1 dT nW(a, T n , {r}), (25) 

a n=0 T„ Jo Jo Jo 

where {r} is a short-hand for the set of times . . . , r„}. The weight is 
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n n 

W(a ! T„,{r}) = (-ir(e-^n e " Tp(Bp "^" l) )(Hll^h)' ( 26 ) 

p— 1 p— 1 

where £" p = (a(p)\D\a(p)) . 

Now, consider an SSE index sequence S n — (pi, . . . ,b n ), containing m indices b p < My, corresponding to m off- 
diagonal and n — m diagonal operators. Removing all the indices b p > My from S n results in a valid sequence T m . We 
use the notation [S n ] for this "projection" of S n onto the corresponding T m ; [S n ] — T m . Since there are no convergence 
issues for a finite lattice model at finite (3, neither for SSE nor for the perturbation expansion, the weights of the two 
formulations must be related according to 

00 f0 f T l l-Tm-l 

V V W(a,S n )= dn dm-- I dT m W(a,T m ,{T}). (27) 

n=m[S„] =Tm J ° J ° J ° 

Hence, for a given sequence of off-diagonal operators, the time integrals of the perturbation expansion correspond 
to a summation over all possible augmentations with diagonal operators in the SSE scheme. The equality ( p7[ ) can 
be explicitly verified in the extreme case where there are no diagonal operators in H, i.e., D — (for example, the 
XY-model in the ^-component basis). In this case [S n ] — S n — T n , and the r-integrals above give (3 n /n\, which is 
exactly the prefactor of the SSE weight (^|). 

For a nonzero D, the dominant SSE strings contain a finite fraction of diagonal operators, and the integrals in the 
interaction representation become nontrivial. In constructing a simulation algorithm based on one of these expansions, 
one hence has to weight the disadvantage of a longer operator string in the SSE scheme against a more complicated 
weight function in the case of the perturbation expansion. Although the perturbation series integrand is formally 
simple, it does not appear to be feasible to carry out the time integrals analytically. It is, however, straight-forward to 
include an importance sampling of the times in the Monte Carlo procedures. It is likely that sampling the perturbation 
expansion will be more efficient than the SSE algorithm in cases where the diagonal term dominates. The average 
length of the perturbation expansion is then significantly shorter than the SSE string. Note, however, that the average 
length of the perturbation expansion has the same scaling ~ (3N as the SSE string length, as will be discussed further 
below. 

In constructing a function A (the estimator) measuring an operator A on the configuration space, symmetries 
of the space should be taken into account in order to reduce the statistical fluctuations. An evident one is the 
translational symmetry of periodic (non-random) lattices. Another one is the periodicity in the imaginary time (or 
SSE propagation) direction, originating from the cyclic property of the trace operation. In the SSE scheme, this is 
manifested by 

W(a,S n ) = W(a(p),S n \p\), V = 0, . . . , n, (28) 

where S n [p] is the index sequence obtained by p times cyclically permuting S n , and a(p) refers to the p times 
propagated state (pT[). One can therefore average the measurements over all p = 0, . . . , n, an example of which is seen 
in Eq. ©. 

The perturbation expansion involves the imaginary times only in the form of differences. The part of the weight 
(Eq) containing the times can be rewritten as 



e Y[e~ Tp(Ep ~ Ep - l] = n e_BpAp ' ( 29 ) 



-(3E 

p=l 



where A p = t p — t p+ i + cr/3 is the time difference between the operators at positions p, p + 1, with a = 0, 1 chosen such 
that A p € [0, (3\, and r n = tq. Shifting all times by an equal amount 8 therefore does not change the weight, provided 
that the shifted times Tj + 6 obey the limits of the time-ordered integration. Hence, S € [— r n ,/3 — t\]. Cyclically 
permuting T n and {r}, and including an appropriate uniform shift S, is also allowed. We define a permutation of 
the times with an implicit shift, such that the last time in the permuted sequence is at its lower bound 0. Denoting 
the p times permuted set {r(p)}, the zero times permuted {r(0)} hence corresponds to a uniform shift 5 = — r n . A 
permutation with an additional shift is denoted {r(p) + 5}. Hence we have 

W(a,T n ,{T}) = W(a(p),T n (p),{T(p) + 5}), <5e[0,A p ]. (30) 

We now derive expressions for some important types of operator expectation values within the perturbation expan- 
sion scheme, and contrast them with those of the SSE formulation. 



G 



First, we consider an equal-time correlation function between two diagonal operators, 



Ci 



(AjAi 



(31) 



In the SSE approach, the expansion of Tr{AjAie /3H } leads to the same sum as in Eq. (||), with each term multiplied 
by the eigenvalue dj [Ola, [01 = (a| A;Ai|a). Using the cyclic property ( ]2q ) then leads to Eq. (|l7|). In the interaction 
representation, Eq. (|30| ) implies that each cyclic permutation p should be weighted with the time interval A p . Since 

E P A p = ft we S et 



C V = ~Q \ \ a i\P\ a j\p\ 



(32) 



\p=i 



This type of expression is, of course, valid for any diagonal operator. 
The SSE expression for the static susceptibility, 

Xa = / driAjWMO)), 
Jo 



(33) 



can be obtained by integrating the time-dependent expectation value ( [15|) over r. Alternatively, one can include a 
source hjAj in the Hamiltonian, and calculate the response function via 



Xij = 



d(Ai 



hi=0 



The result 



iffl 



Xij 



ft 



n(n + 1) 



ft 



\p=o 



\p=o 



(n + 1) 



(34) 



(35) 



p=0 



In the interaction representation, the derivative (Q) applied to (Ai) — J2 P ^p a j\p\/ P gi ves 



X>j 



(36) 



Hence, with both methods, there is a simple exact estimator for the static susceptibility. This is important in view of 
the fact that the discretization can introduce spurious temperature dependences in divergent susceptibilities calculated 
using standard worldline methods, due to a combination of Trotter errors and numerical integration errors.Q 

A second class of observables easily accessible in SSE as well as real-space path integral formulations is one involving 
the operators Ht present in the Hamiltonian. First, consider a single operator (Hf,). In the SSE formalism, the 
estimator is 



(37) 



where N(b) is the total number of indices = b in the sequence S n . This formula can easily be derived by noting 
that the expansion of the numerator Tr{e _/3 - ff i?h} leads to a one-to-one correspondence with a subset of the terms 
in Eq. (||), namely, those for which the last index b n = b. The terms are related by a factor — n//3, which hence is 
the contribution to (Hb) from these partition function configurations. For the terms with b n ^ b the contribution is 
zero. Averaging over all cyclic permutations then gives (37). From this result, Eq. ( fl2"| ) for the SSE internal energy 
estimator follows. 

In view of the relation (|27|) between SSE and the perturbation expansion, we would expect Eq. (^) to be valid 
also for an off-diagonal operator in the interaction representation scheme. Proceeding as in the SSE derivation, there 
is again a one-to-one correspondence between the terms of the expansion of Ti{e~^ H Hb} and a subset of those in 
Z. However, the situation is complicated by the fact that, for a given power n, the terms of Z have one more time 
integration. In order to properly relate the terms to each other, we can formally introduce another integral, 



1 



dr n = 1, 



(38) 
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in the expansion of the numerator. Terms of order n — 1 are then in a one-to-one correspondence with terms of order 
n in the partition function ((2^). The lack of the time-dependent exponential associated with the last operator H b and 
the factor l/r n _i in ( |3g| ) imply a contribution —e Tn ( En ~ En -^ /r n -i, if b n — b. By the cyclic property j30|), this can 
be averaged over the time range A p for each p, giving 



(H b ) 




I b ( P )K(p)\, (39) 



where Ib(p) = 1 if b p — b and Ib(p) — otherwise. Since t p _i — t p + A p _i, the contribution if b p = b is 

K(p)= / dT p = e - A P-i(E p -E p ^) / dx _ (4Q) 



This integral cannot be solved in closed form, except if E p — E p -\ = 0. We expect (|^) and (^) to be equivalent. 
Therefore, the average of ( pp| ) over all times must give 1. In fact, this is the case already for the average over all 
T p G [t p +i, t p _i], which is equivalent to the average over all A p _i in the allowed range [0, A p _i + A p ], with A p _i + A p 
kept constant. In doing this averaging, the integral (Eof) has to be weighted by the relative probability of a given 
A p _i, which according to Eq. ( p9|) is ~ e Ap " 1 ' Bp_£ ' p " 1 ' — The resulting double integral can be solved, with the result 

K{p) = ~ = L ( 41 ) 



Thus, we have shown that (|39j) indeed reduces to the SSE estimator (|3j). This then also implies that the average 
length of the perturbation expansion is given by (n) = /3|(t^)|, which scales as (3N. 

In order to derive an expression for an off-diagonal equal-time correlation function of the type 

F(b u b 2 ) = (H bl H b2 ), (42) 

one can proceed along the same lines as for the single operator considered above. The SSE expression is again formally 
very simple. Each occurrence in S n of a pair of indices hd>2 gives a constant contribution. Denoting by N(b\b2) the 
number of such pairs of adjacent operators, the result isEZI 

F(b 1 ,b 2 ) = {(n-l)N(b 1 b 2 ))/p 2 . (43) 
In the interaction representation formalism, the estimator is 

F(b 1 ,b 2 ) = (^I blb2 (p-l,p)K(p-l,p)J, (44) 

where I blb2 {p — l,p) = 1 if the indices at the adjacent positions p — l,p are b\ and b 2 , and zero otherwise. In order 
to evaluate the contribution K(n — 1, n) from a pair at the last two positions, n — 1 and n, we now insert a double 
integral 

dr n -x \ d,T n = l, (45) 



(Tn-2? 



for a term of power n — 2 in the numerator. Performing the appropriate cyclical permutations and time averages 
analogous to the ones discussed above, the contribution from an arbitrary pair of adjacent operators is 



f p dz f z dv f " dx\e x{E ?- E '-^ 
K(p - 1 p) = Jo az JQ ay Jz ax ^ e ( 46 ) 

/(f" dyJ Dp dxev( E p- E p-i)e x ( E p-i- E P -^ ' 



where 



D p = A p _ 2 + A p _j + A p . (47) 



Calculating the integrals results in 
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K(p -lp) = (E P -i-E p _ 2 )/p 

(B p _ 1 -_B p _ 2 )(e- D P (E f- B P-2)-l) 

where special cases such as E p — E p -2 = should be treated as limiting values. In order to relate this much more 
complicated expression to the simple SSE result of a constant contribution [n— 1)//3 2 [Eq. (|i3|)], we note that a typical 
value of the time interval D p is 3/3/n, and K(p — l,p) ~ l/(/3D p ) for D p small. Hence, a typical value of — l,f>) 
is ~ n/P 2 , i.e. of the same order as the SSE contribution. 
The SSE estimator for the static off-diagonal susceptibility, 

Xtj = / dr(i? b2 (r)F bl (0)), (49) 
Jo 

is given by the remarkably simple formula0 

Xtj = (N^Nih) - Sb^Nih)) If3. (50) 

As this expression only involves counting the numbers of indices b\ and 62 in the sequence, it must [by the configuration 
relation, Eq. (p7|)1 be the correct expression in the interaction representation as well [in contrast, the equal-time 
correlation function discussed above involves pairs of indices, the distributions of which are different in SSE and the 
interaction representation, due to the diagonal operators present in SSE] . We shall not prove this explicitly here. 

The above derivations have clearly shown the close relationships between the SSE configuration space and the 
continuous time path integral. We note that in cases where the expressions differ, they are generally formally simpler 
in the SSE case. In this sense, the SSE propagation dimension is a more natural representation of the quantum 
fluctuations than imaginary time. 

At this stage a reader may wonder why the expansion is dominated by such large powers (n) ~ N(3, independent of 
the size of the perturbation, and why is it not dominated by small powers when the perturbation theory converges for 
an infinite system at T = 0. The answer lies in the fact that the stochastic sampling is done for the partition function, 
for which there is never a convergent expansion as T — > and the size of the system goes to infinity. To clarify the 
situation, let us assume that the off-diagonal part of the Hamiltonian is multiplied by a perturbation parameter A. 
Furthermore, the free energy per unit volume, /, has a convergent expansion 

/(A) = f + fiX + 0{\ 2 ). (51) 

Then the series expansion for the partition function for a system of volume N becomes 

2 = e-Wcce-W>=)]fl ri r, (52) 

n 

with |a„| = (/3-ZV) n |/i|™/(n)!. It is easy to show that a n is maximum for n — /3N\fi\ in agreement with Eq. (|37 

An updating scheme for importance sampling of the perturbation series can now be constructed along the lines 
discussed in the previous section in the context of the SSE method. The differences are only in the weight function, 
which involves a set of times which also is sampled stochastically. As a pedagogical example, in the following section 
we develop the details of an algorithm for the simple case of the anisotropic Heisenberg chain. In Sec. V we extend 
the scheme to include couplings to phonon degrees of freedom (spin-Peierls model). 



IV. ALGORITHM FOR THE HEISENBERG CHAIN 



Here we describe the details of a perturbation series algorithm developed for the anisotropic 5 = 1/2 Heisenberg 
chain. We discuss some properties of the method and use exact diagonalization results for small systems as well as 
known analytical results for the thermodynamic limit to show that very accurate, unbiased results can indeed be 
produced. 



A. Construction of the Algorithm 



The model we consider here is defined by the Hamiltonian 
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N 



H = J^2[SiS? +1 + ±(S+S7 +1 + S? +1 Sr)], (J>0), (53) 

i=l 

with periodic boundary conditions. A controls the anisotropy, with A = 1 corresponding to the isotropic Heisenberg 
point. 

We wish to construct an algorithm in which, as in the SSE scheme, Monte Carlo updates that change the expansion 
order, n, are accomplished by inserting or removing diagonal operators one at a time, and off-diagonal operators are 
inserted or removed pairwise in substitutions with diagonal operators. Since there are only off-diagonal operators in 
the perturbation expansion string, we add constants to the Hamiltonian, and formally consider these as part of the 
perturbation. For the spin chain at hand we define 



H ub = -I (54a) 

i+l°i ' 



% = fc+S+Sr, (54b) 



and write the perturbation as 



V = 



■A- 2 N 



o=l 6=1 



For iV even, only operator strings with an even number of the off-diagonal operators #2,6 contribute to the partition 
function. The weight (^6|) is hence positive definite. The matrix element of an allowed operator string equals one, 
and therefore 

n 

W(a,T n ,{r}) = (A/2) n e- pBo '[[e- T ^ E '- E '- 1 \ (56) 

P =i 

T n is now for convenience defined as a sequence of index pairs 

T n = [01,61], [02,62], ■ ■ ■ , [a n ,b n ], (57) 

with dp <E {1,2} and b p 6 {1, . . . , N} referring to the operator type and lattice bond (nearest-neighbor spin pair), 
respectively. Using the fixed-length scheme developed for the SSE algorithm, we define Hq q — I, and insert L — n of 
these in each string. Taking into account the number of possible insertions then gives 

W(a,T L ,{r}) = (A/2) " ( ^ ~ n)!n! e -^° f[ e -n>(*»-g,-0 (58) 

p=i 

where n is the number of non-[0, 0] elements in Tj,. There are no times associated with the augmentation operators 
Hofi, and the index p in the product therefore refers to the pth non-[0, 0] operator. 

The Monte Carlo sampling is based on updates of the types ( |l9| ) and (|20|). Using [a,b] p as an alternative to the 
notation [a p , b v ], the update changing the power n by ±1 is 

[0,0] p «-> [l,b] p , (59) 

and the simplest update involving off-diagonal operators is 

[l,b] pi ,[l,b] P2 ^ [2,b] pi ,[2,b] P2 . (60) 

These local updates are sufficient for generating all operator strings for an open chain, or a periodic system in the zero 
winding number sector. An update changing the winding number will be discussed further below. For a simulation 
in the canonical ensemble, i.e., with the total magnetization m — ^ i Sf fixed, no further updates of the state |a) are 
required, since ( |60| ) also implies flips of nearest-neighbor spins in the propagated states \a(p)) with p = pi, . . . ,p 2 — 1 
(here and in the following, the periodicity of the sequence is always implied, so that if p\ > P2, the affected states are 
Pi, . . . , L — 1, 0, . . . ,P2 — 1)- In the grand canonical ensemble, global spin flips changing the magnetization of all the 
propagated states also have to be carried out. 

Fig. [l] shows a graphical representation of a configuration generated for an 8-site isotropic Heisenberg chain. This 
type of representation emphasizes that in this simulation scheme the ordered sequence of operators is the central 
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element, and the times formally can be thought of as auxilliary variables associated with the operator positions. In 
the simulation, only one of the states |a(p)) needs to be stored, since all the other ones are uniquely defined given the 
operator sequence and can be generated as needed. 

First, consider the single-operator update (p9[). This can be attempted consecutively at all positions p = 1, . . . , L 



for which a p € {0, 1}. In calculating the Metropolis acceptance probabilities for such updates, the fact that the 
augmentation operators Ho,o are not associated with time integrals has to be taken into account. We now define 
the time difference A p = r(< p) — r(> p), where r(< p), and r(> p) are the times associated with the non-[0,0] 
operator closest to p, with position indices < p and > p, respectively. In a substitution [0, 0] p — » [1, b] p , 6 is chosen at 
random, and a time r p in the range [r(> p), r(< p)] is generated. In the direction [1, 6] p — ► [0, 0] p , the only action is 
to replace [1, b] with [0, 0] and discard the time r(p). One can easily verify that detailed balance with the distribution 



( p8| ) is maintained with the following acceptance probabilities (note that the energy difference E p 
[a p ,b p ] = [1,6]): 



E„ 



p-i 



P([0,0] p ^[l,b] p ) =min 
P([l,6] p -> [0,0] p ) =min 



(A/2)A p iV(n + l) x 
L — n 
L - n + 1 
L(A/2)(A P + A p+1 )iVr 



for 

(61a) 
(61b) 



The pair substitutions ( |60| ) are associated with constraints. In the — ► direction, the first requirement is that 
S£{pi — 1] = — 5£ +1 [pi — 1]. In either direction, an update flips the spins S b [p] and Sfc +1 [p] in the states \a(p)) with 
Pi < p < P2- This implies that there may be no operators [2, 6—1] or [2, 6+1] present between p\ and P2, and hence 
that Sf +1 [p] = —S§ +1 [p] for all p in the range affected by the update. The change of the weight (|58|) in an allowed 
update at 6 then depends on the energy differences E p — P p _i in the local 4-spin substates \S^_ 1: S§, S§ +1 , S b+2 ) p for 
Pi < p < P2 , but only for those p with operators acting on the spins of this substate. _ 

The locality of the constraints and the weight changes allows for a fast updating carried out on subsequencesiH A 
subsequence 6 contains all operators [2, 6] present in Tl, and those [1, 6] operators appearing between antiparallel spins 
S% and S% +1 . Information on the constraints imposed by the presence of the nearest-neighbor operators [2, 6—1] and 
[2, 6+1] is also part of the subsequence. Operators [2, 6—2] and [2, 6 + 2] do not impose constraints, but affect the edge 
spins of the substates \S^_ 1 ,S§,S§ +1 ,S§ +2 ), and therefore also have to be included in the subsequence, so that the 
acceptance probabilities can be calculated. All the 4-spin substates acted upon by the operators of the subsequences 
are also stored for this reason. Clearly, subsequences 6 and 6' can be updated independently of each other if |6' — 6| > 2. 
Since we normally study chains with N a multiple of 4, we simultaneously construct the subsequences for all bonds 
separated by 3 other bonds, and update these one by one. Four such partition-updating cycles are then needed for 
updating all the bonds. 

The subsequence information is arranged as follows: The length of subsequence 6 is denoted Lb, and is the number 
of operators [1, 6], [2, 6], [2, 6—2], and [2, 6+2] present in Tl- These operators are represented by the integers 1 — 4, 
and are stored in lists Ab(l, ■ ■ ■ >£&)• The original positions of these operators in Tl are needed for re- merging the 
updated subsequences into an updated full sequence, and are stored in lists Pf,(l, . . . , Lb). The constraining operators 
[2, 6+1] do not have to be stored. Instead, lists F&(1, . . . , n&) are created, such that Fb(i) = 1 if there are constraining 
operators (one or several) in Tl between positions P,(i) and Pb(i + 1), and Fb(i) = otherwise. The 4-spin substates 
are encoded as single integers (1 — 16), and stored in lists Sb(l, ■ ■ ■ , Lb). 

For updating the subsequences, we use the scheme introduced in Ref. [I?] (other methods are also possible). An 
attempt to carry out a substitution (^0|) in a given subsequence 6 consists of the following steps. A position i\ such 
that Fb(ii) — is first chosen at random. One then searches in the forward direction for the first position j for which 
Ab(J) — Ab(ii) and Fb(j) = 1. This position is the one furthest away from i\ which can be considered, together with 
i\. in a pair substitution. Note that position i = 1 follows i = Lb due to the periodicity, and the search is terminated 
at i = i± — 1 if this position is reached (and then j = i\ — 1). During the search, the positions i of all encountered 
operators Ab(i) = Ab(i\) are stored. One of these, i = 12, is then selected at random, and the pair {Ab(ii), Ab(i2)} is 
replaced by {A' b (ii), A' b (i 2 )} = {2 — Ab(h), 2 — Abfa)} with a probability satisfying detailed balance. 

The total probability of making a certain pair-substitution is the product of the probability P se i e ct[^4(*i)) ^.(^2)] 
of selecting the operators at positions i\ and 12, and the acceptance probability P aC cept[^4(*i)^(*2) - * A' {i^ALli^)]- 
One can show that the selection probabilities with the above procedures are the same in both directions,0 i.e., 
Pscioct[^4(*i), ^1(^2)] = Pseiect[^4'(*i), ^4'(*2)], if the attempt is cancelled with a probability 

p ^^* A ^- NW + k*Y (62) 
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where N(A) and N(A') are the numbers of operators A and A' found between i = i\ and i = j in the search [excluding 
the first operator Af,(ii) = AL If the attempt is not cancelled at this stage, a final acceptance probability is calculated 
on the basis of the weight ( Eq ) , using the stored substates, and the corresponding states modified due to the replaced 
operators (obtained by propagating with the updated subsequence segment). The Metropolis acceptance probability 
is 

P2 

P acccp t[AA -> A' A'] = mm[e-^ E o- E °) [] e-r^^-E'^+E^E^) ^ i] , ( 63 ) 

i=Pl 

where Ei, and E[ are the eigenvalues of D calculated on the substate \S^_ 1 , S^, S^ +1 , S^ +2 )i before and after the 
operator substitution. Typically, a constraint is encountered already a few steps from the starting position i±, and 
therefore the number of operations required per step is rather small (the cancellation probability (|6^) is often 0). 

For each subsequence, a number of updating attempts proportional to the number of operators in the subsequence 
is carried out. The average length of the subsequences, and hence the number of operations needed for its updating, is 
proportional to f3 at low temperatures. After updating all the subsequences belonging to one out of the four partitions, 
the updated operators are inserted in the full index sequence, and changes in the local 4-spin substates are copied 
into the stored full-system state \a). The same procedures are repeated for all four partitions. 

In a periodic system, configurations with a non-zero winding number are possible, and cannot be generated by 
the local updates discussed above. A winding number corresponds to an excess of spin flips in one direction in the 
course of the propagation with the operator string, i.e., a cyclic permutation of same spins in \a{L)) with respect to 
those in |a(0)). The winding numberj-ean be changed by substituting a "half-ring" of off-diagonal operators by the 
complementary half-ring according tct£l 

[2, bi] pi , [2, b 2 ] P2 , • ■ ■ [2, b N / 2 ] PN/2 <-> [2, b' 1 ] pi , [2, b' 2 ] P2 , . . . [2, b' N / 2 ] PN/2 , (64) 

where the bonds 6i, . . . , 6jv/2j b[, . . . , b' N , 2 are a permutation of all the bonds of the periodic lattice. The acceptance 
rate for this type of update decreases rapidly with increasing system size, due to the increasing number of constraints. 
In practice, simulations for systems larger than N « 16 — 20 must be restricted to the sectpc.with zero winding number 
The resulting small error is a boundary effect, and vanishes in the thermodynamic limit r3 

At high temperatures, simulations can be carried out in the grand canonical ensemble, by including updates of 
the total magnetization. This is in principle easily achieved by flipping "straight lines" of spins, Sf [0], . . . , S*[L — 1], 
which is allowed provided that there are no operators [2,i — 1] or [2, z] present in Tl, and the acceptance probability 
then depends only on the neighbor spins S-LjJO] and Sf, 1 [0]. However, the likelihood of an allowed spin flip decreases 
rapidly with decreasing temperature, and in practice simulations for T < J/10 have to be carried out in the canonical 
ensemble. 

Finally, we also perform updates of the times {r} without changes either in the operator string or the states. A 
single time t p can be updated by generating a new time in the allowed range [t p +i, r p _i] (and t\ < f3, r„ > 0), and 
accepting this with a Metropolis acceptance probability calculated from Eq. (p8|): 

P(r P -> r' p ) = min [cxp[(r p - T p )(E p ^ - E p )], l] . (65) 

The typical time separation, and hence the difference t' p — t p , scales as 1/iV, and is typically very small. The acceptance 
rate for these single-time updates is therefore close to 100% in most cases. It is clear that the rate of evolution of 
{t} updated this way [in addition to the generation of a random time when inserting an operator in an update (|59|)], 
will be very slow for large systems. Therefore, we consider simultaneously updating a whole set of times {r Pl , . . . t P2 }, 
using the following scheme. 

A position pi is first chosen at random, and pi is chosen as the smaller of p± + n T and n, with n T a number chosen 
randomly between 1 and some upper bound m T , and m T is adjusted so that a reasonable acceptance rate (« 50%) is 
maintained. If the weight would be independent of the times, the distribution of the times would be uniform within 
the limits of the time ordered integral. If the separation r P2 — t Pi is not too large, the true distribution will be close to 
uniform. We therefore attempt to replace the selected set of times by a randomly generated ordered set {r pi , . . .r p2 }, 
with t' Pi < r pi _i (t 1 < /3 if pi = 1) and r p2 > t P2+ i (r p2 > for pi = n). The Metropolis acceptance probability for 
this multi-time update is 



P{{r P } - {r'» = min 



exp[J2(T p -T p )(E p ^-E p )) ,1 



(66) 



It is clear that the acceptance rate is essentially determined by the time difference t P2 — t Pi , independently of the 
system size. In simulations of large systems, the maximum number of simultaneously updated times, m T + 1, can 
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therefore be as high as 10 3 or higher (in many cases, all times can be updated simultaneously). The importance of 
the multi-time updates will be further discussed below. 

We define a Monte Carlo step (MC step) as a sequence of diagonal updates (^9|) at all positions in Tl, followed 
by off-diagonal pair updates ( p0|) at all bonds. In a grand canonical simulation, a global flip of each spin is also 
attempted, and for simulations of small systems with fluctuating winding number "ring updates" d&j] ) are carried out. 
The number of multi-time updates per MC step is chosen such that, on the average, «50% of the times are changed. 
As already noted, an MC step requires of the order of N(3 operations. 

B. Performance Tests 

We now present some tests of the accuracy of the method. We also briefly address the issues of autocorrelation 
times, and the equilibration of the simulation. We only consider the spin isotropic Heisenberg case [A = 1 in Eq. (p3|)l - 
Detailed numerical finite-temperature results for this model have recently been obtained using the SSE method and 
high-temperature expansions .c3 

In order to verify that the new QMC algorithm indeed produces results free from detectable systematic errors, we 
carried out a long simulation of a 12-site system at inverse temperature f3 = 8, in the grand canonical ensemble and 
with fluctuating winding numbers. These results can be directly compared with exact diagonalization data. 

A useful internal check of the simulation in the isotropic case is the internal energy calculated in two different 
ways; from the diagonal nearest-neighbor correlation function as E\ = 3(S* S* +1 ), according to Eq. (p"7|), and from the 
expectation value of the off-diagonal operators, E% = — (3/2) (#2,1), according to Eq. (p^). For the 12-site system at 
(3 = 8, a simulation consisting of 2 x 10 8 MC steps gave E\ - 0.44372(3) and E 2 = -0.44368(2), where the numbers 
within parentheses indicate the statistical errors. The exact result is E = —0.443697. Hence, both QMC estimates 
are accurate to within relative statistical errors of less than 10~ 4 . 

We also calculated the static structure factor 

S(9) = ^I>~ i0 '~ 0, < 5 i 5 ">' (67) 

3,1 

and the corresponding static susceptibility 

*(?) = jjJ2 e ~ iU ' l)q [ 3 dT{S?(r)S?(0)). (68) 

3,1 

Comparisons with the exact results are presented in Table |. Here the relative accuracy is the highest, sa 10 -4 , close 
to q = 7r/2. The accuracy is lower close to q — tt, due to the strong antiferromagnetic fluctuations present in the 
model. At q — the accuracy is hampered by the low acceptance rate for the spin flips required in the grand canonical 
ensemble. 

The above results clearly demonstrate that the method in principle is very accurate. We now consider a significantly 
larger system, and address some practical issues that arise in realistic simulations. 

For the Heisenberg chain, a number of exact results are known in the thermodynamic limit. For example, Eggert 
et al. recently calculated the temperature dependence of the uniform magnetic susceptibility \ — xil ~ * 0); using the 
thermal Bethe Ansatz.LJ We here compare their infinite-size result with QMC data for a system with 128 spins. The 
simulations consisted of 3 — 5 x 10 7 MC steps for each temperature. At high temperatures {T/J > 0.1) the grand 
canonical ensemble was used, and at lower T the magnetization was kept fixed at m z =0. In the latter case, we 
define the uniform susceptibility as xili) where q\ is the smallest non-zero wavenumber; q\ — 2tt/L. Fig. ^ shows the 
results. The relative statistical errors are « 10 -3 down to about T/J = 0.15, below which the accuracy diminishes 
due to the low acceptance rate for the grand canonical spin flips. More accurate results could have been obtained 
by extrapolating to zero wavenumber on the basis of a few points close to q = 0, instead of just using q = 0. For 
the grand canonical simulations, the agreement with the exact result is clearly very good, indicating no detectable 
effects of finite size for TV = 128 at the temperatures considered. For the canonical ensemble, all results are slightly 
below than the exact curve, suggesting larger finite-size effects when the magnetization is not allowed to fluctuate 
(extrapolating to q — instead of using qi gives a still larger deviation). 

In any Monte Carlo simulation, care has to be taken that the equilibrium distribution is reached before measurements 
are taken. In our algorithm, the cut-off of the perturbation expansion at order L is determined during the equilibration 
part of the simulation, by monitoring the power n and increasing L whenever n exceeds some threshold fraction of 
L (typically 90-95%). We have found that an efficient updating of the times {r} is crucial for achieving a rapid 
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equilibration. In Fig. || we show results for the cut-off versus the simulation time for an N — 128 system at j3 = 8 
for two different ways of updating the times; consecutively updating all the times one- by-one, according to Eq. j65|), 
and collectively updating a whole range of times according to Eq. (fS6|). With single-time updates, the final cut-off 
L = 1064 was reached after approximately 20000 MC steps, whereas the multi-spin updating gave the final L already 
after 288 steps. In this case, all times were updated simultaneously, with an average acceptance rate of s» 70%. 

Measuring the expectation value of the constant operators Hi^ also provides for a good check of the equilibration 
of the sequence length. With J = 1, —2(Hi t b) should equal 1, but before equilibrium is reached a calculation using 
Eq. (^) will deliver a smaller value. In Fig. || we also show the evolution of this measurement with MC time, for 
the two different ways of updating the times. With multi-time updating, the measured values fluctuate around 1 
already after 10 3 steps, whereas with single-time updates the result are too low even after 10 5 steps. In both cases, 
subsequent simulations consisting of 10 7 MC steps gave comparable results for all calculated quantities, with similar 
statistical errors. Hence, in this case the multi-time update appears to be important for reaching equilibrium, but does 
not significantly affect the statistics of the simulation otherwise. For larger system the single-time update becomes 
even less efficient, and it is likely that the multi-time update then will become important also for the statistics of the 
simulation (i.e., the autocorrelation times of measured quantities). 

The efficiency of the updating process in generating independent (uncorrelated) configurations is measured by 
autocorrelation functions. The optimal frequency for measuring expectation values on the generated configurations 
is in principle determined by such functions. For a quantity A, the autocorrelation function is defined as 

(A{ l + t)A^))-(A(i))* 
aA[t) = (AW) (A(i))* ' (69) 

where A(i) denotes the value of the simulation estimator after i MC steps, and the normalization has been chosen such 
that ctu(O) = 1. For large time separations t, one expects a^it) to decay as e~ t ^ A , where £ Q is the autocorrelation 
time. In principle, the long-time dynamics of the simulation is characterized by a single autocorrelation time, i.e., the 
longest {;a- An observable of interest may or may not overlap with the quantity corresponding to this slowest mode. 
Hence, in practice one will deal with different autocorrelation times for different quantities. 

Successive measurements provide independent information only if they are separated by a number of steps sufficiently 
large for the corresponding autocorrelation function to have decayed substantially. In practise, the decay may be much 
faster than e~*^ a for short times, and therefore the relevant separation between measurements may be shorter than 
£ Q . In general one expects the autocorrelations for observables related to the long-distance properties of the model, 
such as S(tt) and x( 7r ) f° r the antiferromagnetic Heisenberg chain, to decay slower than those for essentially local 
quantities such as the energy, or correlation functions away from the wave number of the dominant fluctuations. 
Fig. ^ shows results for the logarithms of the autocorrelation functions of E, S(tt), xO 1 ")) an d x(0) for a 128-site 
system at (3 — 8. For S(ir) and xM) the asymptotic linear decay is the same, as expected, and the autocorrelation 
time is « 40 MC steps. Note, however, tha cism (i) decays slightly faster than a x ( w ) (t) for short times, which can be 
understood on account of S(ir) being an equal time correlation function, whereas x( 7r ) involves also an integration over 
imaginary time. This effect becomes even more pronounced at lower temperatures. For x(0) the autocorrelation time 
is considerably longer at this temperature (w 80 MC steps), due to the low acceptance rate for the global spin flips. 
In contrast, at higher temperatures x(0) exhibits the shortest autocorrelation time. The autocorrelation function 
for the energy decays very rapidly, and it is difficult to see a regime with a linear behavior before statistical noise 
dominates the measurements. One would expect cxe (i) to contain a component of the slowest mode of the simulation, 
but the overlap may be too small to detect (and hence the effect on the statistics is essentially irrelevant). Clearly, 
the relevant time scale for measurements of the energy is much shorter than its asymptotic autocorrelation time. 

In practice, it is not feasible to measure the autocorrelations in every case. As usual, binning the data, so that 
a single bin represents a simulation time much longer than the asymptotic autocorrelation time, ensures proper 
estimates of averages and statistical errors. The only concern then is that the time spent on the measurements should 
not dominate the simulation. Measuring every 10 — 20 MC steps typically only results in a minor overhead. 



V. SPIN-PEIERLS MODEL 



The effects of phonons on spin systems in one dimension are of great current interest, in view of the recent discoveries 
of spin-Peierls transitions in the inorganic compounds GeCuOs (Ref. |2^) and a'— Nav^Os (Ref. |3|). These materials 
consist of weakly coupled spin-1/2 chains and dimerize at T ^ .14 K and 36 K, respectively. A variety of experiments 
have probed their static and dynamic magnetic properties So 

To date, most numerical work on spin-Peierls models in ID has focused on Lanczos exact diagonalization studies 
of the ground state of dimerized chains, and the finite temperature susceptibility has been calculated using complete 
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diagonalization of undimerized chains.E3~E3 The QMC approach discussed in this paper should be an excellent tool for 
detailed non-perturbative studies of spin systems including dynamic, i.e., fully quantum mechanical, phonon degrees 
of freedom. 

Early studies of ID electronic models including phonons were carried out by Hirsch and Fradkin, using the worldlinc 
QMC method with the phonons introduccd-iiia their real-space displacement coordinates, in imaginary time discretized 
with the standard Trotter approximation. EJ We have developed an extension of the Heisenberg algorithm described 
in the previous Section, with phonons treated in the real-space occupation number basis. We believe that this is a 
more efficient way to include the phonons, and furthermore the scheme is exact. The same ideas could also easily be 
implemented within the SSE scheme, but so far we have not done so. The present scheme is likely more efficient than 
SSE at high temperatures, where the average phonon occupation number is high and dominates the energy. At low 
temperatures, SSE may in fact be slightly more efficient. 

We here restrict ourselves to the perhaps simplest type of spin-phonon coupling, namely, we consider a dispcrsionlcss 
harmonic oscillator with frequency ojq at each bond i (Einstein phonons). The exchange interaction between the spins 
at sites i and i + 1 is modulated by the phonon displacement Xi . The Hamiltonian is hence 

N N 

H = y~l( J + Qia:i)Si ■ Sj+i + u>o }^ n». (70) 

i=l i=l 

Here J is the "bare" exchange, a is the spin-phonon coupling, and is the phonon occupation number, given in 
terms of the phonon creation and destruction operators as rn = a\ai. We define the phonon displacement operator as 

Xl = (ot + en)/V2, (71) 

and hence absorb in the coupling a the mass and spring constant of the oscillator. The oscillator potential corre- 
sponding to the definition ( fnj ) is V(x) = ujqx 2 /2. 

We remark that the linear coupling is not completely realistic, since it can lead to negative (ferromagnetic) spin-spin 
couplings if the fluctuations are large. In the low temperature regime this does not occur, and the model should then 
be a good starting pointibr understanding the effects of phonons in real quasi-lD spin systems. 

Augier and PoilblancEj have recently studied a model related to the Hamiltonian (fr0|). They included only the 
phonons with momentum q = tt, which are the ones forming a condensate in the dimerized state, and carried out 
T = Lanczos exact diagonalization in a space with a restricted number of these phonons. They found quantitative 
changes in the behavior from that of a statically dimerized system. Here we will show that the q ^ tt phonons are 
also important at finite temperature, in particular the q — ones which lead to a temperature dependent effective 
spin-spin coupling. 



A. Algorithm 

In a simple modification of the Heisenberg algorithm described in the previous Section, we now have additional 
phonon occupation number operators in the unperturbed Hamiltonian; 

N N 

D = J J2S z b SI +1 +u;o^b- (72) 

6=1 6=1 

We write the perturbation V — J2 a b ^a,b in terms of the following operators: 

(73a) 
(73b) 
(73c) 
(73d) 
(73e) 
(73f) 
(73g) 
(73h) 



Hl,b — 


-Jo/2 






#2,6 = 


(J /2)(5+Vi- 


+" ^b+l^b ) 




^3,6 = 


(a/2V2)(S+S H 


1 + ^b+l^b 


>6 


Hi,b = 


(a/2V2)(S+S^ 


1 + Sb+l^b 


~)a b 


i?5,6 = 


(a/V2)(S* b SI +1 


-!)«! 




^6,6 = 


(a/V2)(S§SI +1 


- |) a 6 




Hr.b — 


{a/2V2)a\ 






Hs,b = 


{a/2V2)a b . 
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The definitions of and He,b contain off-diagonal operators ~ a\ and ~ a^, which are not present in the original 
Hamiltonian. They have been included in order to make the weight function positive definite. They also induce a 
constant shift in the oscillators, which has no effect other than shifting the value of J by an amount proportional to a 2 . 
The operators Hj^ and Hg^ have been included in the Hamiltonian in order to enable straight-forward measurements 
of correlation functions of the displacement operators Xb, which are off-diagonal in th e ch osen r epre sentation and 



hence must be measured using expressions such as Eqs. (|4J) and (|50|) [the prefactor in ( |73g|) and ( 73h ) is chosen for 
convenience, but is in principle arbitrary]. Again, these operators just shift the effective value J and therefore do not 
affect the physics of the model. 

In order to compensate for the shift induced by the added operators, we define the spin-spin coupling in Eqs. ( |73b| ) 
and ((Z2|) with a constant Jo 7^ J, determined such that the new bare coupling J = 1 (i.e., its value in the absence 
of the spin-phonon coupling aJ^jXjSj • Sj+i, which in turn induces an additional shift which is part of the physics 
of the model). For each bond, the added operators equal (3a/4)x&. Thus the shifted oscillator potential is V'(x) = 
u!o(x — 3a/4o;o) 2 /2. The actual bare value of J is then 

J = J + ^-, (74) 

and we therefore choose Jo = 1 — j(a/ J)(a/u)o) and J = 1. If we do not want to deal with negative values of Jo, this 
implies that the maximum coupling constant we can study, for a given w , is a max = 2(Jw /3) 1 / 2 , which is probably 
not a serious restriction for realistic situations. We can relax the restriction by choosing a smaller prefactor in the 
definition of Hjb and -ffs.b, or leaving out these operators altogether (which can be done if no phonon correlation 
functions need to be measured in the simulation) . We note that in principle one can use Jo < and avoid the coupling- 
constant restriction. However, this causes a sign problem due to frustration. The standard way of treating the phonons, 
in the representation where the displacements x% are diagonaljij avoids this by enforcing a positive coupling. This 
corresponds to including non-linear terms in the oscillator potential and/or the coupling, and hence modifies the 
model itself. This clearly makes sense from the point of view of modeling real materials. Including nonlinear terms 
that strictly enforce a positive sign appears to be complicated with our simulation scheme, and therefore the standard 
method is likely preferable in cases where this would be necessary. Note, however, that including "simple" nonlinear 
terms in the oscillator potential should be possible, as long as they do not cause a sign problem. 

It is now a simple matter to include Monte Carlo updates involving the operators H a ,b with a = 3, . . . , 8 in the 
Heisenberg simulation scheme developed in the previous Section. As before, we denote by Hq^ a unit operator, not 
part of the Hamiltonian, introduced in order to fix the length of the index sequence to L. The operator string matrix 
element in the weight function, Eq. ( p6|) , is more complicated than before, since the matrix elements of the phonon 
operators depend on the occupation numbers n&; 

ab\n b ) = \fnb~\n b - 1), 

a\\n b ) = y/n b + l\n b + l). (75) 

Hence, changes in the matrix element of have to be evaluated in pair substitutions involving phonon operators. 
Apart from this, the methods we use for updating the configurations are very similar to the Heisenberg case, and we 
shall therefore merely list the types of substitutions carried out and only very briefly comment on their implementation. 

As before, the update changing the expansion power n is [0,0] <-> [1,6], with acceptance probabilities given by 
Eq. (|6l]) [with Jo in place of A]. The pair substitutions can be divided into two classes. A substitution of the first 
type is only accompanied by changes in the spin states. We use the following ones: 

[1,6], [1,6] ~ [2, 6], [2, 6] (76a) 
[3,6],[4,6]~ [5,6], [6,6] (76b) 
[4,6], [3,6] ~ [6,6], [5,6]. (76c) 

These updates are again carried out with Tl partitioned into one out of four sets of subsequences. The other updates 
cause changes only in the phonon states. These are 

[2,6], [2,6]^ [3,6], [4,6] (77a) 

[3,6], [4,6]^ [4,6], [3,6] (77b) 

[1,6], [3,6] «-> [5,6], [2,6] (77c) 

[1,6], [4, 6]^ [6,6], [2,6] (77d) 

[2, 6], [5, 6]^ [3, 6], [1,6] (77e) 
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[2,6], [6,6] «-> [4,6], [1,6] 
[5,6], [6,6]^ [7,6], [8,6] 



(77f) 
(77g) 



as well as the ones obtained from these by permuting operators both on the left and the right. Since the oscillators are 
independent (i.e., we have not included phonon-phonon couplings), these updates can be carried out independently 
for all bonds, with Tl partitioned into N subsequences. 

For small systems, we also carry out "ring" updates changing the winding number according to Eq. (|64|). Global flips 
of single spins can again be carried out only at relatively high temperatures, thus necessitating the use of a canonical 
ensemble for the spin sector at lower temperatures. For the phonons, the algorithm is automatically grand canonical 
at any temperature, since the phonon number is not conserved. We do, however, also include global fluctuations of 
the number of phonons at each bond. For this purpose, we monitor (during equilibration) the average number of 
phonons per oscillator, (rib), and attempt adding or subtracting between 1 and yl (n&) phonons at a time. 

Since the number of phonons is unbounded (in practice the largest number sampled is limited, but can be very 
large at high temperatures), we cannot test the QMC program against exact diagonalization results with the full 
phonon space included. However, we have compared results with exact diagonalizations for a 4-site chain with the 
number of phonons per oscillator restricted to less than or equal to two. This restricted case already encompasses all 
the elements of the simulation algorithm, and hence we can conclude that our program generates results exact within 
statistical errors. 



We shall here only present some selected illustrative results for the spin-Pcicrls model, deferring to a separate 
publication a detailed discussion of— the physics of the model, and its relevance in interpreting various experimental 
results for quasi-lD spin systems.E-3 The main purpose here is to give an impression about the kind of detailed 
information that can be extracted using the new QMC algorithm. We do, however, also point out some important 
model features of experimental relevance that can be inferred already from the results obtained here. 

We note that in the adiabatic limit ujq — > 0, the model should dimerize for any ao The situation for finite frequency 
is still not completely settled. For the case of noninteracting spinless fermions, which is equivalent to the XY spin 
chain, Hirsch and Fradkin found numerical evidence that a critical phonon coupling is required for a dimerization to 
occur. We here consider a phonon frequency ujq — J/10, and a coupling a — J/4. Our low-temperature results show 
that the system is dimerized at T — for these parameters. Based on their exact diagonalization study including the 
q = 7r phonons only, Augier and PoilblandfJ suggested a coupling a/ J w 0.38 (after adjusting for a difference y/2 in 
definitions) and a bare frequency uoq/J = 0.3 for describing a'— NaV2 05. Hence, the case considered here is closer 
to the adiabatic limit, and the spin-phonon coupling is slightly weaker. This is still a regime that can be expected 
to be of relevance for real materials. As we will see, our bare coupling J is in fact renormalized to a slightly higher 
value due to q — phonons, so that our effective phonon frequency and coupling are somewhat lower in units of the 
effective spin-spin coupling. 

We present simulation results for a wide range of finite temperatures, including low enough T to give the ground 
state for all practical purposes. Systems with up to 128 sites are considered, and a simulation typically consisted of 
1 — 2 x 10 7 MC steps. At high temperatures (T/J > 0.125), we used the grand canonical ensemble for the spins, and 
at lower T we had to restrict ourselves to the canonical ensemble with m z — 0, due to the low acceptance rate for 
global spin flips. For the relatively large lattices considered, the effects of this restriction are minor. As already noted, 
the simulation algorithm is automatically grand canonical for the phonons, and we have found no practical problems 
with the ability of the phonon numbers to fluctuate effectively. 

The spin-phonon coupling causes an average uniform displacement (a;,-) > of the oscillators, due to the energy 
lowering realized by increasing the average coupling, J — > J + a(xi) — J c fr, balanced by the increased potential energy 
of the oscillator. Hence, the spin-Peierls system is characterized by a temperature dependent effective spin-spin 
coupling J c ff(r). Fig. |5| shows our QMC results for the effective coupling, along with its RMS fluctuation, defined as 



As already noted, the relevance of the model to real materials is likely limited to the regime where J e g rarely fluctuates 
to negative values. The results of Fig. ||| indicate that, for the model parameter considered here, this is the case for 
T/J < 1, where J ofr > a(J cS ). 

At low temperatures, the lonq-wavelength spin susceptibility shows the behavior expected for a dimerized spin 
chain. Fig. ^ shows results for x(l) f° r 1 < t/2. The behavior changes rapidly between inverse temperatures (3 = 16 



B. Results 




(78) 
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and 64, from having only a very weak q-dependence in the q — > regime in the former case, to a rapid suppression 
in the latter case. This clearly shows the presence of a spin gap. There is only a minor change going to (3 = 128, 
indicating that the system is essentially in its ground state at this inverse temperature. In order to accurately extract 
the uniform susceptibility x(Q — > 0), it is necessary to extrapolate to q = using 2 — 4 low-q points. This procedure 
is of sufficient accuracy in the interesting intermediate temperature regime where the susceptibility drops rapidly. 

The uniform susceptibility is one of the most important experimentally measurable quantities. It is often used 
to estimate the Heisenberg spin-spin coupling. In comparing results for the spin-Peierls model with the Heisenberg 
chain, it is natural to measure the temperature in units of the effective coupling at T = 0, which we denote J c ff(0). 
We here have J c ff(0) « 1.273 [see Fig. |[. Susceptibilities calculated in the simulation, where the bare exchange J = 1 
is used to set the temperature scale, then also have to be adjusted by this factor since the definition ([38|) contains 
the inverse temperature. Hence, unless stated otherwise, wc now define (3 = J c g(0)/T, and give the susceptibility in 
units of l/J o ff(0). |— - 

In Fig. |?] we compare the uniform susceptibility of the spin-Peierls model with the exact Heisenberg result eB There 
is a significant shift in both the peak position and the amplitude. Both can be understood on the basis of the reduced 
antiferromagnetism due to the fluctuations induced by the phonons. There are no signs of a suppression of the 
susceptibility for T / J c g(0) > 0.1, below which there is a sudden, very rapid drop. Above this drop, we find that the 
shape of the temperature dependence of \ can be quite well fit to the result for the Heisenberg chain without phonons, 
but with an exchange Jm < J e ff(0), as also indicated in the figure. The over-all magnitude of the susceptibility is then 
found to be lower than for a Heisenberg chain with this exchange J = Jgf In terms of an effective Landee g-factor 
(i.e., the value of the ^-factor one would deduce under the incorrect assumption that the system is described by a 
Heisenberg chain), which is 2 for an ideal 5 = 1/2 Heisenberg chain, this implies g$x < 2. For the parameters used 
here, we find Jg t ~ 0.82 J and ga t « 1.86. This is interesting, since experimentally it is [Sometimes found that the 
(/-factor extracted from the susceptibility fit is less than 2, whereas one would in fact expeclca a g- factor slightly larger 
than 2. For example, Eggert was able to fit well experimental susceptibility data for Sr2Cu03 with the Heisenberg 
x(T),E3 but the g-factor corresponding to the fit is as low as 1.6, and this seeming discrepancy has not been explained. 
Based on our results presented here, we propose that the apparent g-factor reduction may be at least partially due 
to phonons. Additional calculations for various phonon frequencies and couplings are underway to further investigate 
this possibility. Note that the susceptibility of the spin-Peierls model deviates significantly from the Heisenberg x(T) 
also if one considers a temperature-dependent energy scale J e s(T) instead of the fixed scale J c s(0) used above, as 
shown in the inset of Fig. [7[ The shift in the peak position remains, and there is a regime where the phonons cause 
a significant susceptibility enhancement. This shows that the effects of the phonons cannot be simply captured by 
a mapping to a Heisenberg chain with a temperature-dependent exchange equal to the average coupling J c ff(T) of 
the modeLi.e., the fluctuations in J c ff(T) must be taken into account. These issues will be discussed in more detail 
elsewhere. c3 

The effects of the phonons on the antiferromagnetism can be seen directly in the staggered susceptibility, graphed 
in Fig. ^ For the Heisenberg chainp^(7r) = Dln 1 ^ 2 (A/T)/T at low temperatures, where the constants have been 
estimated as D w 0.32 and A ps 6.c3 For the gapped spin-Peierls model we should have xW — > constant. We 
find here that there is a maximum in xM approximately at the same temperature at which the rapid drop in the 
uniform susceptibility is seen [the staggered structure factor S(ir) exhibits a similar maximum]. This is behavior is 
distinctively different from that of a gapped system with fixed values of the exchange couplings, such as a statically 
dimerized chain, for which x( 7r ) simply saturates at a value set by the gap. In the presence of dynamic phonons 
the staggered correlations among the instantaneous spin-spin couplings only build up gradually as T is lowered, and 
therefore the antiferromagnetic spin correlations can initially grow stronger than what they will eventually be in the 
statically dimerized T — state. Experimentally, this may have implications for NMR experiments, which probe 
the low- frequency dynamic spin-spin correlations. The NMR rates for models including phonons should be accessible 
using the same numerical techniques (QMC and maximum-entropy analytic continuation) that have recently been 
used for the standard Heisenberg chain.E3 : E3 

The phonon correlations are also interesting, and can give some indirect information on the dynamics as well. We 
define the static phonon structure factor and susceptibility according to 

3,1 

= il2 e ' i{j ~ l)q f dT ^) - <*»(^(0) - <*»>• (79b) 

iV 3,1 J ° 

We expect both to develop peaks at q = n at low temperatures, signaling the dimerization instability. Fig. |^ shows 
results at several temperatures. At high temperatures, both S x (q) and Xx(q) are almost independent of q and increase 
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with increasing temperature. This reflects the behavior of (almost) independent harmonic oscillators. At lower 
temperatures, the expected peaks at q = tt develop. In the high and intermediate temperature regime (T > loq) 
graphed in Fig. |[ Tx x {q) — S x (q) within statistical errors. Hence, at these temperatures the phonon dynamics 
is for all wavelengths dominated by the first Matsubara frequency and is thus essentially classical. Note that the 
phonon susceptibilities have considerably smaller statistical errors than the structure factors. This is contrary to 
the situation for the spin structure factor and susceptibility, which are calculated using the diagonal S* operators. 
This counterintuitive feature of the QMC method can be traced to the simple form of the off-diagonal susceptibility 
estimator, Eq. (|50|). 

The static phonon structure factor and susceptibility can be related to the dynamic real-frequency phonon correla- 
tion function (or spectral function) A(q,w) through sum rules. We define the spectral function as 

/oo 
dte' iuJt (xi(t) Xj {0)). (80) 

From the Lehmann representation one can derive the following sum rules in the standard way:@ 

S x (q) = - (l + e-^)A(q,w) (81a) 

I" JO 

**(«) = - / {l-e-^)-A{q,w). (81b) 

At T = 0, these sum rules can be used to obtain an upper bound for the lowest phonon excitation of momentum q, 
according to 

w min (q) < 2S x (q)/ Xx (q). (82) 

In the case of a single sharp phonon mode, 2S X {q)/Xx (q) is the exact excitation energy-. This bound has recently 
been used to extract spin and charge velocities of ID electronic models from QMC data. El Here we demonstrate its 
use for the phonon spectrum. Fig. |l^ shows results for a 128-site system at a low temperature, T — J/128. For long 
wavelengths, the calculated bound is at the bare phonon frequency luq/ J = 0.1 within statistical errors, indicating a 
very weak effective coupling of these phonons to the spin system. For q — > tt, there is a clear reduction, as expected due 
to the softening of the q = tt mode in a system that spontaneously dimerizes. In an infinite system, the T — > bound 
2S X (tt)/Xx (it) — > since a static dimerization implies the classical relation Xx{^) — (3S x (tt). For a finite system, there 
is always a gap (decreasing with increasing size) to the lowest phonon excitation and both S x (q) and \x (q) therefore 
saturate at finite values for all q. For q ^ tt, the results shown in Fig. [l(] are saturated, but x(7r) (but not S x (tt)) still 
grows with decreasing T and the actual T = bound is therefore considerably lower than that in the figure. 

We now calculate the size of the Peierls distortion. In the dimerized state the average displacement alternates 
between even and odd sites; (xi) — (x) ± S. Hence, the T — staggered phonon structure factor is, for large system 
sizes N, given by 

S x (ir) =a 2 5 2 N. (83) 

Fig. [n] shows S x (ir) versus the inverse temperature J/T for several system sizes. As expected, the saturation takes 
place at a temperature which decreases with increasing system size. The results for N = 32 and N = 128 obey 
Eq. (|||) quite accurately, the ratio of the saturated values being w 4. For N — 8 the value is considerably lower, 
indicating much larger fluctuations in this small system. Using the N = 128 result with Eq. ( p3[ ) gives a dimerization 
aS w 0.13, i.e., the effective alternating average exchange is Ji = J c ff(0)[l ± Aj], with Ay w 0.10. 

The main purpose here has been to illustrate how various physical quantities are accessible with the new QMC 
algorithm. We have therefore not discussed how ouuesults compare to, e.g., mean-field theory and previous numerical 
calculations including only the q — it phonon modeo These important issues will be addressed in a future publication. 



VI. SUMMARY AND DISCUSSION 



In this paper we have introduced a new quantum Monte Carlo algorithm based on the standard perturbation expan- 
sion in the interaction representation. This starting point was first suggested by Prokof'ev et aZJla Our implementation 
of the sampling of the series is different-. and is essentially an adaptation of procedures previously developed for the 
Stochastic Series Expansion algorithm.EJ We have shown that the SSE sum and the continuous imaginary time path 
integral are in fact very strongly related to each other. 
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As in the SSE scheme, the central element of the new method is the ordered sequence of operators. In the interaction 
representation formulation developed here, time enters in the form of auxilliary variables associated with the sequence 
positions. The operator updates leading to changes in the particle propagation paths (kink-antikink creation and 
annihilation, in the language of Prokof'ev et alxJi) are carried out via substitutions of pairs of off-diagonal and 
constant operators, with the time field held fixed. An efficient procedure for collective updating of whole segments of 
the time field was introduced. We also derived expressions for several types of important operator expectation values, 
and compared these with the corresponding expressions previously obtained within the SSE scheme. 

Not having carried out comparisons with the procedures suggested by Prokof'ev et ai, we do not know which of the 
two sampling schemes is more efficient — most likely this is model dependent. We believe that the method discussed 
here is more efficient than the SSE approach in cases where the diagonal part of the Hamiltonian dominates the 
internal energy. However, the weight function is simpler in the SSE case, and again it is not easy to make a general 
statement of exactly when the advantages of new method become significant. 

We note that the loop-cluster algorithm, invented bypEvertz et aZ.,c3 has recently been used with considerable 
success in studies of various spin-1/2 Heisenberg systems,E3 in particular in the exact formulation developed by Beard 
and Wiese.lii, Variants of the loop algorithm have also been suggested for the ID Hubbard modelj£3 and for S > 1/2 
spin chains.c^l However, physics results for these cases have not yet been presented. Although its performance for the 
Heisenberg model can be considered-spectacular (however, for some quantities more accurate results have actually 
been obtained with the SSE methodca), it is not clear how the loop algorithm will fare with more complicated models, 
such as spin-phonon models. It is known to break down completely in some cases.EHI In contrast, the stochastic series 
expansion algorithm and the perturbation series schemes are in principle completely general and are practically useful 
for a wide range of models for which the sign problem can be avoided. 

As a demonstration of the power of the new method, we have implemented it for studies of a spin-Peierls model. 
We considered oscillators associated with the bonds, coupled to the_spins via a linear modulation of the exchange. 
Unlike earlier studies of ID electronic models coupled to phonons,EJ we used the occupation number basis also for 
the phonons. The new QMC algorithm is ideally suited for this type of models, since the bare phonon part of the 
Hamiltonian is diagonal. The initial results presented here for the spin-Peierls chain indicate that reliable results can 
be obtained with modest computer resources (the simulations were run on mid-range workstations). The method 
should be very useful for resolving issues related to the effects of dynamic phonons on the physics of the recently 
discovered. inorganic spin-Peierls compounds, as well as other quantum spin systems. Work along these lines is in 
progress.ti] ID itinerant electronic models of the Hubbard and t-J types including phonons can also be studied using 
the procedures developed here. We also believe that studies of spin-phonon models in higher dimensions are feasible 
with our new method. 

The results presented here for the spin-Pcicrls model already indicate some important consequences of dynamic 
phonons at finite temperature. The type of phonons considered here naturally lead to a temperature dependent 
effective spin-spin coupling. We found that the uniform magnetic susceptibility still has a shape rather similar to 
that of the Heisenberg chain in a sizable regime close to the susceptibility maximum often used to extract the size 
of the exchange coupling from experimental data. However, both the value of J and the (/-factor extracted from 
a fit are reduced relative to a Heisenberg chain with a coupling equal to the average coupling of the spin-phonon 
model. We propose thajUjthis dynamic effect may at least partially be the reason for the reduced g- factor found, in, 
some quasi-lD systemsa Furthermore, these results cast some doubts on the validity of detailed extractions£§EZI 
of the nearest-neighbor and next-nearest-neighbor couplings in GeCu03 from fits of exact diagonalization data of 
frustrated Heisenberg chains to susceptibility measurements. It is likely that the couplings extracted from such fits 
do not directly correspond to the true spin-spin couplings of the system, but are influenced by the temperature 
dependence of the couplings as well as their fluctuations, as discussed above. Interchain couplings likely also have 
some non-negligible effects on the susceptibility. Unfortunately, the QMC approach introduced here does not allow 
for studies of frustrated systepSr-Xat least not at very low temperatures), due to sign problems. Since a'-NaV2C>5 is 
not expected to be frustrated,c3E-3 we believe that detailed experimental studies of this material in combination with 
finite-T QMC studies will be of key importance in clarifying the microscopic physics of the spin-Peierls materials. 
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g/vr S(g) (QMC) S(g) (exact) x (q) (QMC) x{q) (exact) 

0.008581(24) 0.008590 0.06864(20) 0.068720 
1/6 0.048110(5) 0.048120 0.11822(3) 0.118266 
1/3 0.102865(8) 0.102875 0.14422(3) 0.144236 
1/2 0.169906(13) 0.169917 0.19746(5) 0.197478 
2/3 0.263758(18) 0.263768 0.32669(8) 0.326648 
5/6 0.43371(5) 0.433742 0.8027(3) 0.803118 

1 0.95472(18) 0.954566 3.9897(13) 3.989397 



TABLE I. Simulation results for the static structure factor and the static susceptibility of a 12-site Heisenberg chain at 
/9 = 8, compared with the exact results. The numbers within parentheses indicate the statistical errors (defined as one 



standard deviation of the averages), i.e. 0.123(45) stands for 0.123 ± 0.045. 
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FIG. 1. A configuration generated for an 8-site periodic system with A = 1 at /3 = 4. A row represents a spin state \a(p)), 
with p — to p — L from top to bottom. Solid and open circles indicate up and down spins, respectively. Dashed and solid 
bars represent operators [1, b] and [2, 6], and their associated times are graphed to the right. 

FIG. 2. Upper panel: The uniform magnetic susceptibility (in units of 1/J) of a 128-site Heisenberg chain calculated using 
the new QMC method (circles with error bars) compared with the exact result for the infinite system (curve). Lower panel: 
The relative deviation (xqmc — Xcxact)/Xcxact of the QMC data from the exact result. The solid and open circles are for grand 
canonical and canonical simulations, respectively. 

FIG. 3. Upper panel: The truncation L vs. the number of MC steps performed during the initial parts of two simulations 
of a 128 site system at /3 = 8. In one case (solid curve), all the times were simultaneously updated, whereas in the other case 
(dashed curve) the times were consecutively updated one-by-one. 

FIG. 4. The autocorrelation function for E (solid squares), x(0) (open squares), x( 7r ) (solid circles), and S(n) (open circles), 
obtained in a simulation of a 128-site system at inverse temperature /3 = 8. 

FIG. 5. The temperature dependence of the effective spin-spin coupling (upper panel), and the size of its RMS fluctuation 
(lower panel). 

FIG. 6. The spin susceptibility (in units of 1/J) vs. wave-number in the long-wavelength regime for a system with N = 128 
at inverse temperatures (3 = 16 (solid circles), /3 = 32 (open circles), /3 = 64 (solid squares), and (3 — 128 (open squares). The 
statistical errors are at most the size of the symbols. 

FIG. 7. The uniform susceptibility [in units of 1/J e fj(0)] vs. temperature [in units of J ef f(0)], compared with the exact 
Heisenberg result (solid curve). For the Heisenberg chain J e g = J = 1. The dashed curve is the Heisenberg susceptibility 
for J = 0.86 x Jefi(O) and a p-factor ~ 1.82. The inset shows the QMC data graphed on the temperature scale set by the 
temperature dependent coupling J e g(T) [\ here also contains this factor, i.e., it is given in units of l/J c g(T)], compared with 
X(T) for the Heisenberg chain. 

FIG. 8. The staggered susceptibility [in units of 1/J e ff(0)] vs. temperature (solid circles), compared with results for the 
Heisenberg chain (open circles). The curve indicates the asymptotic divergent form for the Heisenberg chain. 

FIG. 9. The static phonon structure factor vs. wave number at T/J — 1.0 (open circles), T/J = 0.5 (solid circles), T/J — 0.25 
(open squares), and T/J = 0.125 (solid squares). The corresponding susceptibilities, multiplied by T/J, are indicated by the 
(in some cases barely visible) solid curves. The statistical errors of the susceptibilities are considerably smaller than those of 
the structure factors, and are not indicated for sake of clarity. The dashed lines indicate the g-independent structure factors of 
independent harmonic oscillators at the corresponding temperatures. 

FIG. 10. Upper bound (in units of the bare exchange J) for the lowest phonon excitation vs. momentum for a 128-site system 
at T — J/64. The dashed line indicates the bare, momentum independent phonon frequency luq/J — 0.1. 

FIG. 11. The staggered phonon structure factor vs. inverse temperature for systems of size N = 8 (solid circles), N = 32 
(open circles), and N = 128 (solid squares). The statistical errors are smaller than the symbols. 
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Fig. 1 . A. W. Sandvik, R. R. P. Singh, and D. K. Campbell 




Fig. 2. A. W. Sandvik, R. R. P. Singh, and D. K. Campbell 
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